We look at the difference of symptom rates (percent column) between 2019 and 2020, and consider if there is any significantly long or anomalous time interval showing consecutive lower rates in 2020.

We construct the hypothesis testing:

\[H_0: \text{There is no significant time interval showing lower symptoms reporting rate in 2020 compared with 2019;} \] \[H_1: \text{There exists significant time interval that has low symptom reporting rates in 2020 compared with. } \]

We use the one-dimensional SAT scan statistics to construct the significance. In other words, the time interval is considered to be significant, if it has low probability to be reached by random appearence (probability <0.05).

Testing for Case 1.

Record the difference of reporting rates for case 1 between 2019 and 2020 as ts1.

rm(list = ls())
setwd('C:/Users/andre/OneDrive/Desktop/prometh')
library(readr)
text2020 <- read_csv("C:/Users/andre/OneDrive/Desktop/prometh/dailytext2020.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_double(),
##   days = col_date(format = ""),
##   sent = col_double(),
##   responders = col_double(),
##   ili1 = col_double(),
##   ili2 = col_double(),
##   ili2a = col_double(),
##   ili3 = col_double(),
##   resp.prop = col_double(),
##   ili1percent = col_double(),
##   ili1percentsma = col_double(),
##   ili2percent = col_double(),
##   ili2percentsma = col_double(),
##   ili3percent = col_double(),
##   ili3percentsma = col_double(),
##   ili2apercent = col_double(),
##   ili2apercentsma = col_double()
## )
text2019 <- read_csv("C:/Users/andre/OneDrive/Desktop/prometh/dailytext2019.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_double(),
##   days1 = col_date(format = ""),
##   sent = col_double(),
##   responders = col_double(),
##   ili1 = col_double(),
##   ili2 = col_double(),
##   ili2a = col_double(),
##   ili3 = col_double(),
##   resp.prop = col_double(),
##   ili1percent = col_double(),
##   ili1percentsma = col_double(),
##   ili2percent = col_double(),
##   ili2percentsma = col_double(),
##   ili3percent = col_double(),
##   ili3percentsma = col_double(),
##   ili2apercent = col_double(),
##   ili2apercentsma = col_double(),
##   days2 = col_date(format = "")
## )
#text2019 = read.csv(file = 'dailytext2019.csv')
#text2020 = read.csv(file = 'dailytext2020.csv')

text2020<-text2020[order(as.Date(text2020$days, format="%Y-%m-%d")),]
text2019<-text2019[order(as.Date(text2019$days2, format="%Y-%m-%d")),]

#they both end in 05/21
#but t2019 starts in 1/30
#t2020 starts in 1/27
#cut off t2020 to also start in 1/30
#remove the 1st 3 rows
#adjusting for leap year, 2019 shifts a day ealier, so remove last row in 2020
text2020<-text2020[-c(1,2,3,nrow(text2020)),]

text2019$day = as.Date(text2019$days2,format ="%Y-%m-%d")
text2020$day = as.Date(text2020$days,format ="%Y-%m-%d")
ts2019 = text2019$ili1percent#[order(text2019$day)]
ts2020 = text2020$ili1percent#[order(text2020$day)]

#ts2020_n = ts2020[-c(1,2,3,34)] # remove days 2020 differ from 2019
ts1 = ts2019-ts2020 #difference

Use the median as cutoff. The longest time interval with difference > median is 13 which is located after spring break (starting date: 03-28).

library(ggplot2)
library(scales)
## 
## Attaching package: 'scales'
## The following object is masked from 'package:readr':
## 
##     col_factor
ts1_binary = ts1>mean(ts1)
#hut<-qplot(ts(ts1),ylab = "rate2019-rate2020")
hut<-data.frame(ts=ts1,day=text2019$day)

#h<-median(ts1)
rle.ts = rle(ts1_binary)
max(rle.ts$lengths[rle.ts$values==TRUE])
## [1] 18
#max(rle.ts$lengths[rle.ts$values==F])
pos = which(rle.ts$lengths==18)
#neg= which(rle.ts$lengths==11)
## starting date of this time period
start.pos = cumsum(rle.ts$lengths[1:(pos-1)])[pos-1]
v = c(start.pos+1,start.pos+18) #period 
text2019$day[v] #period
## [1] "2020-04-01" "2020-04-18"
this_theme<-theme(
        plot.background=element_rect(fill="#ffffff"),
        panel.background = element_rect(fill = "#FDF1E6"),
        panel.grid.major = element_line(colour = "lightgrey"),
        plot.title = element_text(hjust = .5, face = "bold", size=20),
        plot.subtitle = element_text(hjust = .5, face="bold", size=17),
        panel.border = element_blank(),
        axis.line.x = element_line(size = 0.5,linetype = "solid"),
        axis.line.y = element_line(size = 0.5,linetype = "solid"),
        legend.text=element_text(size=15,face="bold"),
        axis.text.x = element_text(size = 10,angle=80, vjust=0.6,face = "bold"),
        axis.text.y = element_text(size = 10,face="bold"),
        axis.title.x = element_text(size=15,face = "bold"),
        axis.title.y = element_text(size=15,face="bold"),
        legend.spacing.x = unit(1.0, 'cm'))

p <- ggplot(hut,aes(y=ts1*100,x=day)) +
  geom_line() +
  scale_x_date(breaks=date_breaks("7 days"),date_labels = "%m-%d")+
  ylab("Rate difference \n (2019 - 2020)")+
  xlab ("Date")+
  labs(subtitle = "Fever and either sore throat or cough")+
   geom_hline(yintercept = mean(ts1)*100,colour="red", linetype="solid")+
  geom_vline(xintercept = as.Date(text2019$day[start.pos]+1),colour="blue", linetype="dashed")+
  geom_vline(xintercept = as.Date(text2019$day[start.pos]+18),colour="blue", linetype="dashed")+
    ggtitle("") +
    this_theme

p

Perform permutation test: the longest time interval is significant according to permutation test (p value=0.0036). Hence, the 2020 has a significate time period with lower reporting rates compared with 2019, which is located after spring break.

set.seed(10)
n=10000
rs_max = rep(0)
for(i in 1:n){
  pre1 = ts1_binary[1:112]
  rs = pre1[sample(1:112, size=112, replace=FALSE)]
  rle.rs = rle(rs)
  rs_max[i] = max(rle.rs$lengths[rle.rs$values==TRUE])
}
sum(rs_max>=18)/n
## [1] 1e-04

Check that whether there is any significant time interval before spring break. The longest time interval before spring break is 5, which is not significant (p value 0.8747).

rle.ts$lengths[rle.ts$values==TRUE]
##  [1]  1  2  2  2  5  1  2  4  1  3  2 18  1  4  2  2  1  2  5
sum(rs_max>=5)/n
## [1] 0.9497

Case 2

ts2019_2 = text2019$ili2percent#[order(text2019$day)]
ts2020_2 = text2020$ili2percent#[order(text2020$day)]

#ts2020_n = ts2020[-c(1,2,3,34)] # remove days 2020 differ from 2019
ts2 = ts2019_2-ts2020_2 #difference

Use the median as cutoff. The longest time interval with difference > median is 13 which is located after spring break (starting date: 03-28).

ts2_binary = ts2>median(ts2)
#hut<-qplot(ts(ts1),ylab = "rate2019-rate2020")
hut2<-data.frame(ts=ts2,day=text2019$day)

h<-median(ts2)
rle.ts2 = rle(ts2_binary)
max(rle.ts2$lengths[rle.ts2$values==TRUE])
## [1] 12
#max(rle.ts$lengths[rle.ts$values==F])
pos = which(rle.ts2$lengths==12)
#neg= which(rle.ts$lengths==11)
## starting date of this time period
start.pos2 = cumsum(rle.ts2$lengths[1:(pos-1)])[pos-1]
v2 = c(start.pos2+1,start.pos2+12) #period 
text2019$day[v2] #period
## [1] "2020-03-31" "2020-04-11"
p2 <- ggplot(hut,aes(y=ts2*100,x=day)) +
  geom_line() +
  scale_x_date(breaks=date_breaks("7 days"),date_labels = "%m-%d")+
  ylab("Rate difference \n (2019 - 2020)")+
  xlab ("Date")+
  labs(subtitle = "Sum of symptoms' score")+
   geom_hline(yintercept = mean(ts2)*100,colour="red", linetype="solid")+
  geom_vline(xintercept = as.Date(text2019$day[start.pos]+1),colour="blue", linetype="dashed")+
  geom_vline(xintercept = as.Date(text2019$day[start.pos]+12),colour="blue", linetype="dashed")+
    ggtitle("") +
    this_theme


p2

Perform permutation test: the longest time interval is significant according to permutation test (p value=0.0036). Hence, the 2020 has a significate time period with lower reporting rates compared with 2019, which is located after spring break.

set.seed(10)
n=10000
rs_max2 = rep(0)
for(i in 1:n){
  pre2 = ts2_binary[1:112]
  rs2 = pre2[sample(1:112, size=112, replace=FALSE)]
  rle.rs2 = rle(rs2)
  rs_max2[i] = max(rle.rs2$lengths[rle.rs2$values==TRUE])
}
sum(rs_max2>=12)/n
## [1] 0.0069

Check that whether there is any significant time interval before spring break. The longest time interval before spring break is 5, which is not significant (p value 0.8747).

rle.ts2$lengths[rle.ts$values==TRUE]
##  [1]  1  1  3  7  3  5  1  1  4 12  7  2  1  2  2  1  2  1 NA
sum(rs_max>=7)/n
## [1] 0.4638